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We describe an algorithm which directly determines the quintessence potential from observational 
data, without using an equation of state parametrisation. The strategy is to numerically determine 
observational quantities as a function of the expansion coefficients of the quintessence potential, 
which are then constrained using a likelihood approach. We further impose a model selection 
criterion, the Bayesian Information Criterion, to determine the appropriate level of the potential 
expansion. In addition to the potential parameters, the present-day quintessence field velocity is 
kept as a free parameter. Our investigation contains unusual model types, including a scalar field 
moving on a flat potential, or in an uphill direction, and is general enough to permit oscillating 
quintessence field models. We apply our method to the 'gold' Type la supernovae sample of Riess 
et al. (2004), confirming the pure cosmological constant model as the best description of current 
supernovae luminosity-redshift data. Our method is optimal for extracting quintessence parameters 
from future data. 



PACS numbers: 98.80.-k 

I. INTRODUCTION 

Quintessence, a scalar field slowly rolling on its poten- 
tial, remains one of the most attractive possibilities for 
explaining the observed acceleration of the Universe (for 
reviews, see Ref. Q). A key goal for future observational 
programs is to seek definitive evidence for variation in the 
dark energy density with redshift, which would exclude 
a cosmological constant. In that event, one would then 
seek an optimal determination of dark energy properties 
in the hope of relating them to fundamental physics. 

In this paper we assume from the outset that single- 
field quintessence remains a viable description of obser- 
vational data, i.e. that it has successfully passed tests 
against other dark energy paradigms. Our aim is then to 
obtain optimal constraints on the quintessence potential. 
We do this by passing directly between the quintessence 
potential and the observable quantities, focusing in this 
paper on the luminosity-redshift relation of type la su- 
pernovae (SNIa). We parametrise the potential, and then 
constrain those parameters, along with global properties 
of the Universe, via a likelihood analysis. Additionally, 
we use model selection criteria in order to select the pre- 
ferred level of parametrisation of the potential. 

Although a variant on the general scheme of recon- 
struction, our approach is distinct from those already in 
the literature [2| in that we do not rely on a parametri- 
sation of the dark energy equation of state, which then 
must be related to the dark energy potential via rela- 
tions which may be approximate (see Guo et al. 3] for 
relations in some particular cases). The work closest in 
spirit to our own is that of Simon et al. , who consider 
an extremely general action and expand the quintessence 
potential in Chebyshev polynomials (in the redshift range 
of available data) . They relate the expansion coefficients 
to the redshift evolution of the matter density and Hub- 
ble parameter. Those quantities are then extracted from 
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observations and processed into constraints on the po- 
tential, those constraints however being as a function of 
redshift rather than scalar field value. Their treatment 
is roughly analogous to the inflationary reconstruction 
method whereby observables such as the spectral index 
and tensor amplitude are obtained from data, and then 
related to the inflationary potential via the slow-roll ap- 
proximation 1^. Our present paper is analogous to the 
direct inflaton potential reconstruction method proposed 
by Grivell and Liddle ^] , where the observed power spec- 
tra are predicted numerically directly from the inflation 
potential. 

II. FORMALISM 

Our set-up is relatively straightforward. We assume 
that the quintessence field (j) has a potential V{(j)), which 
we expand as a power series 

1/(0) = Uo + Vi(j) + U20' + V3^^ + -- - , (1) 

where the field is measured in Planck units and (with- 
out loss of generality) we take cj) to be presently zero. 
Note that when we fit these parameters, even in the 
case of "complete" data we do not necessarily obtain the 
MacLaurin expansion of the true potential to the same 
order as this is generally not the best polynomial fit over 
an interval (in the least-squares or minimax sense). We 
choose not to use a Chebyshev series as in Ref. Q for 
the following reasons. Firstly, when fitting, the coeffi- 
cients of a Chebyshev series would just be linear com- 
binations of the coefficients of a monomial basis polyno- 
mial of the same order, so the fits are the same. Fur- 
thermore, Chebyshev polynomials would depend on the 
range (j) takes which in turn depends on the model param- 
eters. Lastly, because we do not fit the potential expan- 
sion directly, but rather through a function depending on 
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an integral of the potential, expanding the potential in 
orthogonal polynomials will not guarantee uncorrelated 
coefficients. 

The scalar field obeys the equation 



dV 



(2) 



with the Hubble parameter given by the Friedmann equa- 
tion 



SttG 



(3) 



Here pm is the matter density and = (p^ /2 + V{(j)) 
the quintessence density. We assume spatial flatness 
throughout, though the generalisation to the non-flat 
case would be straightforward. Since then + = 1 
we have the initial condition 



•0 = ±J2 [(1 - rtndPc.O - Vo] 



(4) 



We allow (/>o to take either sign but results are symmetric 
under simultaneous reversal of its sign and of odd-order 
expansion coefficients. 

In this article we focus on SNIa data, and hence the 
observational quantity we need to predict is the lumi- 
nosity distance as a function of redshift. The luminosity 
distance is given by 



dL{z;e) 



Ho 



(5) 



where 

l + z 



1/2 



(6) 

is the Hubble-constant-free luminosity distance (for the 
low redshifts we consider there is no contribution from 
radiation), and 

F{z-@)^i( (l + w^(z';0))dln(l-Kz'). (7) 
Jo 

is the parameter vector describing the model, and 
_ 4>V2 - V{<i>) 
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P0 02/2 + F(0) 
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is the quintessence equation of state. The apparent mag- 
nitude m{z; 0) of a type la supernova can be expressed 
as 



771(2; 0) ^ M + 5\og 



10 



Mpc 



25, 



(9) 



where M is the absolute magnitude of SNIa (suppos- 
ing they are standard candles). The distance modulus 
0) is defined as 



/i(z;0) = m(z;0) -M. 



(10) 



In the following, we should in principle keep M as a free 
parameter. To this end, we define our "observational" 
quantity to be 



(11) 



where rrii is a measurement of TO(2:i;0true) (with 0truc 
the projection of the parameters of the "true" model of 
the Universe onto our model and its parameter space). 
Supernovae observations measure m, but typically report 
a distance modulus 



H* EE m, - M* 



(12) 



where M* is some estimate of the absolute magnitude. 
The relevant quantity for fits is thus 

IJi - p{zi; 0) = fj,* +T] - 5logiQ'Di^{zi; 0) , (13) 

where 



77 = 5 logio {Ho Mpc) + AM - 25 



(14) 



and AM = M* — M. Because of the perfect degeneracy 
between M and logj^Q Hq , and the fact that the equa- 
tions are otherwise independent of these parameters, our 
effective D-dimensional parameter vector is 



0-(77,0o,K),...,"^^D-3). 



(15) 



Recall that (po and Vq determine the matter density 
through Eq. igj. 

To end this section, we note that our formalism in- 
cludes some possibilities which are not commonly con- 
sidered. Even if the potential is truncated as a constant, 
the present field velocity remains a free parameter and so 
the scalar field can move on this flat potential.^ To regain 
the cosmological constant case we must make the addi- 
tional assumption that this velocity is zero. Further, the 
field may be rolling uphill; this may seem unlikely but 
is valid phenomenologically and might occur in models 
where the field has recently passed beyond a minimum. 
Our analysis can also generate models where the scalar 
field has undergone one or more oscillations about a min- 
imum in the recent past. 



III. DATA ANALYSIS 

A. Likelihood analysis 

We carry out a likelihood analysis of the models in 
comparison to the observational data from Riess et al. j7| . 
We use the 157 SNIa of the 'gold' sample. 



^ This situation is equivalent to having a pure cosmological con- 
stant plus a stiff fluid with w = 1. 
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With a prior distribution P{@), the posterior proba- 
bihty of the parameters 0, given the data set, is given 
according to Bayes' theorem by 

P(0|data) = -|/:(data|0)P(0) = ^e-^'^®'^/^P{@) , 

(16) 

where 

^2(0)^^(M,^i^f^ (17) 

i=l i 

is summed over ah N data points, and 
Z — J £(data|0)P(0)d0 is a normahsation 
constant, irrelevant for parameter fitting (but as it is the 
Bayesian evidence, highly relevant for model selection 
as discussed below). Here and tXi are the observed 
distance moduli and their standard deviations, Zi the 
redshift of the observed supernova and n{zi;@) the 
distance modulus predicted for the redshift Zi by our 
model with parameters 0. 

To estimate parameters we wish to find P(0|data) ex- 
plicitly as a function of 0. This is in general non-trivial, 
and the standard approach is to explore the parameter 
space in some way and keep a histogram characterizing 
P(0|data). We choose to explore the parameter space 
using an MCMC approach 8, 9, El [ill. MCMC calcula- 
tions are generally preferable over grid methods as they 
scale approximately linearly with the dimension of the 
problem, rather than exponentially. 

Our MCMC algorithm is the following, and makes 
use of relatively standard step optimisation and conver- 
gence/mixing testing. 

1. The starting points for the Markov chains are chosen 
to be close to the expected high-likelihood region with 
some random spread, checking that they satisfy the pri- 
ors. 

2. Starting with an initial best guess for the covariance 
matrix of the underlying distribution, we optimise the 
step sizes of the Gaussian trial distribution with the it- 
eration rule 10 

Ct. = (2.4Vi?)C,_i , (18) 

where C^j is the i^^ estimate of the covariance matrix 
of the trial distribution, D is the number of parameters 
and Ci_i the covariance matrix of the {i — 1)**^ chain 
produced (with Cq our initial best guess) . We use chains 
of 10 000 elements for the optimisation process, and con- 
tinue updating the trial distribution until there is no sig- 
nificant increase in the sampling efficiency (assessed by 
comparing the eigenvalues of the covariance matrices). 
Between each iteration, the parameter space is rotated to 
the eigenspace of the new covariance matrix, to maximise 
the efficiency in exploring the shape of the likelihood dis- 
tribution. 

3. The full production run is started. A set of m chains 
with n elements each is generated, and only these are 



used for the final analysis. We generate well separated 
starting points as before for each of the chains. The 
chains are tested for convergence and mixing using the 
Gelman-Rubin test [TH [l^ . which compares the vari- 
ances within a chain to the variances between chains, 
which in the asymptotic limit should give a Gelman- 
Rubin ratio R = I. We require R < 1.05 for each param- 
eter. A consistently high and non-convergent Gelman- 
Rubin ratio is indicative of a very loosely constrained 
parameter. 

In the above, all calculations of covariances and means 
are done by first dropping an initial burn-in section 
from the chain. We define the burn-in section following 
Tegmark et al. 13] as the elements in the chain from 
the beginning up to the first element to have a like- 
lihood value above the median likelihood value of the 
whole chain. The chains were analysed using a slightly 
modified version of GetDist provided with CosmoMC [l0| 
(again with burn-in sections excluded). 

We impose two important constraints on the behaviour 
of the cosmology within the redshift range < z < 2, 
and hence as priors on the parameters. Firstly, the to- 
tal energy density of the universe must remain positive 
at all times to exclude collapsing epochs, and secondly 
we need to avoid models where the kinetic energy would 
dominate at early epochs (such domination is permitted 
by the SNIa data alone, but is inconsistent with other 
data as discussed later). We limit the kinetic contribu- 
tion to Okin < 0.5 for z >1 — see Section FlV CI Most 
marginalised posterior likelihoods are fairly insensitive 
to the particular choice of upper limit. However, some 
marginalised posteriors involving (j)Q do change signifi- 
cantly. This will be discussed further in the Results sec- 
tion. Additionally we naturally impose < ilm < 1. No 
other priors (e.g. on Hq) were found to be necessary to 
obtain acceptable cosmologies. 

B. Model selection 

The order of the power series of the potential can be 
freely chosen, and the results obtained will obviously de- 
pend on the order to which it is taken, with parameters 
becoming less and less constrained as the order increases. 
In addition to a determination of the best-fitting parame- 
ters within a given model, one therefore needs to compare 
the different models (i.e. expansions to different orders) 
in order to determine which is the preferred fit to the 
data. 

Since models with more parameters will always lead to 
an improved best-fit model, one must use model selection 
statistics [ll[ll[ll. These set up a tension between the 
number of model parameters and the goodness of fit. In 
the context of Bayesian inference the best such statistic 
is the Bayesian evidence 0, ; for an application to 
SNIa data see Ref. ^3 ■ The evidence is however difficult 
to calculate, and in this paper we use a sirnpler statistic, 
the Bayesian Information Criterion (BIG) 0,0], which 
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FIG. 1: One and two-dimensional likelihood distributions for 
D = 2. Solid lines are marginalised ID likelihoods and dotted 
lines mean ID likelihoods. Solid 2D contours represent 68% 
and 95% regions of the marginalised distribution, and shading 
reflects the mean distribution. 

gives a crude approximation to the evidence. The BIC is 
given by 

BIC = -2\nC„,^^ + D\nN, (19) 

where /^max is the likelihood of the best-fitting parame- 
ters for that model, D the number of model parameters, 
and TV the number of datapoints used in the fit. Models 
are ranked with the lowest value of the BIC indicating 
the preferred model. A difference of 2 for the BIC is re- 
garded as positive evidence, and of 6 or more as strong 
evidence, against the model with the larger value [13. H^ . 

It is worth mentioning that although we specifically 
consider a quintessence scenario, a model selection result 
favouring more than one potential parameter would indi- 
cate a dynamical dark energy component more generally, 
since for every choice of {H{z), p^aiz)} there exists a cor- 
responding quintessence potential, by virtue of Picard's 
existence theorem for ODE's (demonstrated explicitly in 
e.g. Ref. )• 

IV. RESULTS 

The maximum likelihood value and parameters were 
estimated using the approach described in the preceding 
Section. We investigated the cosmological constant case 
and then cases of one, two and three potential parameters 
(i.e. polynomial orders zero, one, two). Since solving the 
necessary ODEs is not computationally intensive, we can 



generate very long chains. For each scenario, 10 chains 
each containing 1 000 000 elements were obtained. 

A. Parameter estimation 

1. Cosmological constant (D = 2) 

As a check, we investigate the case of a cosmological 
constant. Our parameter vector is 

&^iv,Vo), (20) 

since (pQ — 0. Indeed we obtain the well-known results 
for SNIa data, as seen in Fig. ^Thc constraint on the 
matter density is shown along with that of other models 
in Fig.H 

2. Constant potential with kinetic energy (D = 3) 

Allowing a non-zero kinetic contribution on a constant 
potential means our parameter vector is 

& = {t^,^o,Vo). (21) 

This does, unsurprisingly, improve the fit to data relative 
to the pure cosmological constant. However, looking at 
the likelihood distributions in Fig. O we clearly see that 
00 = is not excluded at a statistically-significant level. 
The bimodality in the 0o distributions is due to the model 
depending only on 0q. 

Since a non-zero kinetic contribution is preferred by 
the data, we also require a higher Vq than in the cosmo- 
logical constant case. A simple way to see this should 
be the case is by considering the effective quintessence 
equation of state: with a kinetic contribution which in- 
creases with redshift, the potential term must be larger 
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FIG. 3; As Fig.[T]for D ^ Z. 
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FIG. 5: As Fig.IUfor D = 4. 
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(a) Posterior 
distribution. 
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(b) Prior distribution. 
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FIG. 4: D = ?>. Posterior and prior distributions for Vo and 

00- 



than for the cosmological constant case to maintain the 
same effective equation of state at high redshift. The 
corresponding shift and spread in U,^ is shown in Fig. |21 
The hmits on (f)Q (and also the other parameters) are 
dependent on our choice of prior on Hkin, but the above 
conclusions remain even in the (unrealistic) case of no 
prior. Hence, the choice of prior on l^kin can be effectively 
regarded as a choice of upper limit on |(/)o|. The cut-off 
of the prior on fikin is illustrated in Fig. ^ 



3. Linear potential (0 = 4) 

For a linear potential, the parameter vector under con- 
sideration is 
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(a) Posterior 
distribution. 
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(b) Prior distribution. 
Shaded region is 
forbidden. 



(22) 



FIG. 6: D = 4. Posterior and prior distributions for Vi and 

00- 



The likelihood distributions (see Fig. show a strong 
degeneracy between (/)o and Vi, which is the main new 
feature compared to I? = 3. This is because a particular 
value of (f> at some earlier redshift can be attained by 
adjusting either (pQ or Vi . Consequently, the best-fit value 
for (f)Q is less than for D = 3, but with a non-zero Vi. 

We also note the bimodality in the 0o ~ Vi distribution. 
This reflects the symmetry under simultaneous change of 
sign of (j)Q and odd-order expansion cocfhcicnts mentioned 
in the introduction. Just as for the case D = 3, the 
prior on fJkin cuts off the likelihood distribution in a high- 
likelihood region (see Fig. (B)). 

An investigation of the linear potential with a different 
emphasis can be found in Ref. 
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The likelihood distribution is an even function of (po- No confidence hmits can be 
given for 0o, as the prior on f2kin cuts in the high-likelihood region and (/)o = is not 
excluded at the 68% level. The upper limit on \(po\ thus corresponds to the maximum 
allowed value according to our choice of prior on f2kin, as discussed in the text, the 
quoted number being the best-fit. 
^ The likelihood distribution is symmetric under simultaneous change of sign of ^0 and 
odd-order potential expansion coefficients. No confidence limits can be given for (j)o or 
Vi for the same reason as above. 

Because of the difficulty in obtaining a convergent /well-mixed sampling, as discussed 
in the text, we choose not to give confidence limits for D — 5. 

TABLE I: Best-fit model parameters and BIC values. Note that these are the best-fit parameter values and confidence limits 
derived from the full D-dimensional likelihood distribution, not the marginalised distributions. 



4- Quadratic potential (D ~ 5) 

The quadratic potential model has the parameter vec- 
tor 

0-(r;,,^o,"^^o, 1^1,^2). (23) 

In this case we find that the third potential parameter, 
V2 , is unconstrained by the data, characterised by a large 
and oscillating Gelman-Rubin ratio (around 1.3-1.9) for 
that parameter. Because of that we do not show the 
likelihoods, since the marginalised distributions will not 
have the correct weights. 

To explore this situation further, we ran additional 
chains using V2 = arctan(V2) as our parameter instead 
of V2 (this corresponds to a change in prior on V2, since 
dV2/dV2 is a function of V2). This choice was motivated 
by the expectation that a cosmological constant, which is 
achieved in either of the limits V2 = or V2 = 00, is very 
likely to be a good fit, and hence allows us to explore 
the infinite range in V2 that might be needed. With this 
choice we get convergent and low Gelman-Rubin ratios, 
and arctan(V2) essentially unconstrained. At V2 Pc.o 
we find a small peak in likelihood, but because the distri- 
bution remains high and nearly flat outside this peak it is 
not possible to constrain the parameter without further 
data. 



B. Model comparison 

We compare the different models using the BIC, which 
uses the maximum likelihood achievable by each model. 
The parameters, likelihoods and BIC values are given in 
Table D 



1. Cosmological constant (D ~ 2) 

The cosmological constant forms the base model for 
our model comparison, and as is well known provides a 
good fit to the data. Indeed, the BIC ranks it as preferred 
to our other models. 



2. Constant potential with kinetic energy (D = 3^ 

As mentioned in the context of parameter fitting, 
00 = is not excluded at a statistically-significant level. 
Including this parameter does allow a somewhat better 
fit to the data, but the BIC penalises its extra parame- 
ter leaving the pure cosmological constant model as the 
preferred description of present data. 

The best-fit D — i cosmology, shown in Fig. [7| ex- 
hibits very strong evolution in from kinetic to poten- 
tial domination over the redshift range of available data. 
This is not entirely surprising: even a tiny kinetic contri- 
bution at present will correspond to a much higher kinetic 
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FIG. 7: Dynamical evolution in the best-fit cosmologies. The graphs of the potential show the range of 4> out to 2: = 2. 



component at earlier times, as only Hubble friction will 
work to decrease (/) in this model. 



3. Linear potential (D = 4) 

As clearly seen, {0o = 0, Vi = 0} is well within the 
preferred region, and so the cosmological constant model 
is embedded within the allowed parameter space of our 
extended model. Accordingly, the model comparison of 
Table U prefers the pure cosmological constant model, 
with the BIC difference arguing quite strongly against 
the inclusion of the two extra parameters. 

The best-fit cosmology (Fig. Cj) is practically indistin- 
guishable from the best-fit for D = 3. This arises from 
the strong degeneracy between (po and Vi; it turns out 
that present data are not discriminating in the orthogo- 
nal direction. 

Within the context of model comparison, a curious 
point to note is that a field rolling on a linear potential is 
quite strongly disfavoured as compared to a field rolling 
on a constant potential (Table QJ. This is because the 
inclusion of a potential slope hardly improves the best- 
fit at all, while costing an extra parameter. This may 
seem quite artificial, but is the conclusion of our phe- 
nomenological approach. One should note however that 



the BIC comparison addresses only how well the differ- 
ent models fit the data at hand; when interpretting as 
a model probability one should bear in mind that this 
conclusion could be overturned if one felt that the prior 
model probabilities were quite different. 



4- Quadratic potential (D = 5J 

As Table 12 shows, the quadratic model is strongly dis- 
favoured by the BIC. Note that although we have not 
necessarily obtained a convergent distribution, we can as- 
sess the model as the BIC only depends on the maximum 
likelihood value. However, with such a broad distribution 
the BIC is expected to be a poor approximation to the 
Bayesian evidence. 

A significant feature of the quadratic potential model 
is that the best-fit has fim ~ 0.05, see Fig.[71 This is of 
course in stark contradiction with many other datasets. 
The evolution of is quite different from that for D — 3 
and D — 4, with the field starting high up on the poten- 
tial, rolling past the minimum and reaching the turning 
point by the present time. However the preferred param- 
eter region includes models with much more reasonable 
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FIG. 8: Marginalised posteriors for D = 3 and D = 4 de- 
pending on the choice of prior on fikin- 



C. Choice of prior on fikin 

As mentioned above, we limit the kinetic contribution 
at z > 1. This is necessary because the SNIa data alone 
favourthe kinetic energy to dominate at z = 1, and by 
inference to be completely dominant at higher redshifts. 
This is in contradiction to almost any other cosmological 
dataset (for instance, the mere existence of high-redshift 
galaxies, and of the cosmic microwave background) , and 
so external priors are necessary to keep us in the physical 
regime. 

The precise choice of this upper limit on Okin is some- 
what arbitrary, and does have a non-negligible impact 
on the posterior distribution. In terms of marginalised 
distributions, the distributions involving 00 show the 
strongest dependence, see Fig. |51 Looking at parame- 
ter estimation, as mentioned above the main impact is 
broader confidence limits. For our model selection anal- 
ysis, conclusions remain unchanged as the cosmological 



constant model is still the preferred model even without 
a prior on fikin for all cases. 



V. CONCLUSIONS 

We have described and implemented a reconstruction 
scheme for quintessence potentials from data, using an 
MCMC likelihood approach, which we applied to SNIa 
data. Additionally, we describe the application of model 
selection using the Bayesian Information Criterion, and 
point out the generality of any positive evidence found 
for dynamical dark energy in this approach. 

As might be expected, the data provides positive evi- 
dence in favour of a cosmological constant in our setup, 
based on model selection by the BIC. A similar conclusion 
was previo usly reached by Saini et al. |l7j and by Bas- 
sett et al. amongst a set of models parametrised by 
different equation of state evolution. Some of our distri- 
butions do however exhibit broad non-Gaussian regions, 
which merits a more detailed model selection investiga- 
tion using the full Bayesian evidence (since the BIC is 
only a reasonable approximation for sharply-peaked uni- 
modal distributions). 

The low-dimensional models are quite well constrained 
by current data. However, we find that if one allows 
power series of order two or higher, the parameters (in- 
cluding the linear one) become unconstrained by current 
SNIa data given our very loose priors. This agrees with 
what is suggested by the analysis by Maor and Brustein 
in the context of distinguishing potential classes. 

We plan to extend this work in several directions. As 
mentioned above, it would be interesting to extend the 
model selection to use the full Bayesian evidence, though 
it will then be essential to consider the issue of prior pa- 
rameter ranges, which the BIC sidesteps. Generalisation 
to include further cosmological datasets is desirable, es- 
pecially CMB anisotropics, though that will require the 
potential expansion to be valid over a much wider range 
of redshifts. Finally, it would be interesting to explore 
whether in the context of tracking models it might be 
possible to eliminate the parameter (po, which ought to 
be determined via an early tracking regime. 
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